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Abstract 

The approach to curve implicitization through Sylvester and Bezout resultant ma¬ 
trices and bivariate interpolation in the usual power basis is extended to the case 
of Bernstein-Bezoutian matrices constructed when the polynomials are given in the 
Bernstein basis. The coefficients of the implicit equation are also computed in the 
bivariate tensor-product Bernstein basis, and their computation involves the bidi¬ 
agonal factorization of the inverses of certain totally positive matrices. 
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1 Introduction 


When studying rational plane algebraic curves, there are two standard ways 
of representation, the implicit equations and the parametric equations. The 
intersection of two curves is more easily computed when we have the implicit 
equation of one curve and the parametric equations of the other, and hence it 
is very important to be able to change from one representation to another. 

We will concentrate on the implicitization problem, that is to say, on finding 
an implicit representation starting from a given rational parametrization of 
the curve. 

In [12] we have presented an approach to the implicitization problem based 
on interpolation using the usual power basis for the corresponding space of 
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bivariate polynomials. However, very recent work [2], related to polynomi¬ 
als expressed in the Bernstein basis has showed the importance of evaluating 
resultants from Bernstein basis resultant matrices directly, avoiding a basis 
transformation. In this sense, in [2] it is indicated that for numerical com¬ 
putations involving polynomials in Bernstein form it is essential to consider 
algorithms which express all intermediate results using this form only. 

Although those papers study univariate polynomials, it must be observed that 
the construction of the resultant matrices can be extended to the case in which 
the entries of the resultant matrix are polynomials. It must also be taken into 
account that the Bernstein basis has also important advantages in the context 
of tracing implicit algebraic curves [13]. 

So our aim is to use bivariate interpolation for obtaining in the bivariate 
tensor-product Bernstein basis the implicit equation of a plane algebraic curve 
given by its parametric equations in Bernstein form (which is the usual situ¬ 
ation in the case of Bezier curves). Although we present all the details with 
an example in exact rational arithmetic, it must be taken into account that 
the process can also be carried out in (high) finite precision arithmetic. In 
that situation some important results of numerical linear algebra we use will 
have a major importance. More precisely the total positivity of certain matri¬ 
ces will be an important issue, as it happens in several instances of computer 
aided geometric design (see, for example, the recent work [18] and references 
therein). 


The rest of the paper is organized as follows. In Section 2 several basic re¬ 
sults will be presented. In Section 3 we introduce the interpolation algorithm 
for computing the implicit equation as a factor of the determinant of the re¬ 
sultant matrix, while in Section 4 we consider some results related to total 
positivity which will be relevant for the solving of the linear system associated 
with the interpolation problem. Finally, in Section 5 we briefly examine the 
computational complexity of the whole algorithm. 


2 Preliminaries 

Let P(t) = (x(t),y(t)) be a proper parametrization of a rational plane al¬ 
gebraic curve C, where x(t) = and y(t) = and gcd{u\,v\) = 

gcd(u 2 ,v 2 ) = 1. A parametrization P(t) = ( x(t),y(t )) of a curve C is said 
to be proper if every point on C except a finite number of exceptional points 
is generated by exactly one value of the parameter t. It is well known that 
every rational curve has a proper parametrization, so we can assume that the 
parametrization is proper. Several recent results on the properness of curve 
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parametrizations can be seen in [17]. 

In connection with the implicitization problem, the following theorem [17] 
holds: 


Theorem 1. Let P = ( x[t ) = y[t ) = be a proper rational parametri- 

zation of an irreducible curve C, with gcd{u\,V\) = gcd(u 2 ,v 2 ) = 1. Then the 
polynomial defining C is Res t (ui(t) — xvi(t), u 2 (t) — yv 2 (t)) (the resultant with 
respect to t of the polynomials Ui(t) — xv\(t) and u 2 {t) — yv 2 {t)). 


Our aim is to compute the implicit equation F(x,y ) = 0 of the curve C by 
means of polynomial interpolation, which taking into account Theorem 1 is 
equivalent to compute Res t {u\(t) — xv\(t),u 2 (t) — yv 2 (t)). 

First of all, we remark that the concept of interpolation space will be essential. 
The following result, also in [17], shows which is in our case the most suitable 
interpolation space: 


Theorem 2. Let P = ( x[t ) = y(t ) = yfpjO be a proper rational parametri- 

zation of the irreducible curve C defined by F(x,y), and let gcd(ui,Vi) = 
gcd(u 2 ,v 2 ) = 1. Then deg y (F ) = max{deg t (ui), deg t (vi)} and deg x (F ) = 
= max{deg t (u 2 ), deg t (v 2 )}. 


Theorem 2 tells us that the polynomial F(x, y) defining the implicit equa¬ 
tion of the curve C belongs to the polynomial space II nim (x, y), where n = 
max{deg t (u 2 ), deg t (v 2 )} and m = max{deg t (ui),deg t (vi)}. The dimension of 
II n ,m( x ,y) is (n + 1 )(m + 1), and a basis is given by {x l yi\i = 0, • • • ,n; j = 
0, • • • , m}. Moreover deg x (F(x,y )) = n and deg y (F(x,y)) = m, and therefore 
there is no interpolation space Ll riS (x, y) with r < n or s < m such that F(x, y) 
belongs to H rjS (x,y). 

Let us note that these theorems refer to the degree of polynomials in the 
power basis, so since now we will be using the Bernstein basis some care will 
be needed. For the sake of clarity we will illustrate all our results with a small 
example. Let 


{/3< 4| («)./3! 4 b).d 4) («).d 4 b).d 4) w} 

be the (univariate) Bernstein basis of the space of polynomials of degree less 
than or equal to 4, where the Bernstein polynomials are defined as follows, 

# ) W = (")(!- ; = o,...,n. 
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and let us consider the algebraic curve given by the parametric equation 


, _ 4/3q \t) + 4(3[ \t) + 3/^2 \t) + 3/?3 \t) + 7(3j \t) 
Po 4 \t) + Pi 4 \t) + p2 4 \t) + $\t) + 3/?i 4) (t) 

y(t) — 2/?q \t) + 3 p[ \t) + 3/?2 ^(f) + 3/?3 ^(£) + 4/?4 •’(t). 


If we call p(f) = Ui(t) — xv\{t) and g(£) = u^t) — yv2(t), their coefficients in 
the Bernstein basis are given by 

Po = 4 — x, pi = 4 — x, P 2 = 3 — x, pz = 3 — x, p 4 = 7 — 3x, 

and 

<?o = 2 - y, q 1 = 3 - y, q 2 = 3 - y, q 3 = 3 «• y, q A = 4 - y, 


However, let us observe that 

p(t) = 4 — x — 6t 2 + 8t 3 + (—2x + l)f 4 
(a polynomial of degree 4 in t), while 


q(t) = 2 — y + At — 6t 2 + 4t 3 , 
a polynomial of degree 3 in t. 

Therefore, the polynomial defining the implicit equation will be a polynomial 
belonging to the space U n m (x, y) with n — 3 and m = 4. We will use for that 
space the tensor-product bivariate Bernstein basis given by 

{B-j ,m) } = {Pi n \x)$ m \y),i = 0,... ,n; j = 0,... ,m}. 


Finally we will recall, following [2], the algorithm for constructing Bernstein- 
Bezout matrix of p(t) and q(t). Although in [2] the coefficients of the polynomi¬ 
als are always numbers, in our application we will construct the symbolic (i.e. 
with the entries being polynomials in x, y) Bernstein-Bezout matrix of p(t) 
and q(t) which we denote by BS. For the reader’s convenience, we present the 
algorithm written in Maple language: 


for i from 1 to n do 

BS [i, 1] : = (n/i)*(p[i]*q[0]-p[0]*q[i]); 
od; 

for j from 1 to n-1 do 
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BS[n,j+1]:= (n/(n-j))*(p[n]*q[j]-p[j]*q [n]) 
od; 

for j from 1 to n-1 do 
for i from 1 to n-1 do 

BS [i, j+1] : = (n~2/(i*(n-j)))*(p[i] *q[j]-p [j] *q[i]) 
+((j*(n-i))/(i*(n-j)))*BS[i+1,j]; 
od; 
od; 


Let us observe that if m = n, the resultant is the determinant of the Bernstein- 
Bezout matrix, while -as a consequence of the corresponding result for the 
Bezout resultant [16]- if m > n, that determinant is equal to the resultant 
multiplied by the factor (p m (t )) m ~ n , where p m (t) is the leading coefficient of 
p{t) in the power basis. 

So in our example, the determinant oi BS will be the implicit equation we 
are looking for multiplied by the factor (— 2x + 1), since the degree of p is 
4 and the degree of q is 3 and the coefficient of f 4 in p is (—2x + 1). In the 
following section we will show how to compute the coefficients in the bivariate 
tensor-product Bernstein basis of the implicit equation (which will be a scalar 
multiple of the resultant computed by using the approach of [12], where the 
equation is obtained in the usual power basis). 


3 The interpolation process 


Since the expansion of the symbolic determinant is very time and space con¬ 
suming, our aim is to compute the polynomial defining the implicit equation 
by means of Lagrange bivariate interpolation, but using the Bernstein basis 
instead of the power basis. A good introduction to the theory of interpolation 
can be seen in [5]. 

If we consider the interpolation nodes (.) (i = 0 ,■■■ ,n;j = 0, ••• , m) 
and the interpolation space II n m (x,y), the interpolation problem is stated as 
follows: 

Given (n + l)(m + 1) values 


fij e K (i = 0, • • • , n; j = 0, • • • , m) 
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(the interpolation data), find a polynomial 


E(x,y) = Y, c ijPi n \ x )Pj m \y ) e n n ,m(x,y) 
(i,j) el 


(where / is the index set / = {(i,j)\i — 0, • • • , n; j — 0, • • • , m}) such that 


F{ x i,y j ) = f ij V (■ i,j)ei . 

If we consider for the interpolation space H ntm (x,y) the basis 

i = 0,..., n; j = 0,..., m} = 

{/3} n) (z)/?j m) (y ), * = 0, • • •, n; j = 0,..., m} = 

1^00 ? -£*01 5 * * * 5 ^0m 5 -£*10 5 -£*11 7 * * * ? -£*lm 7 * * * 7 

p>(n,m) T?(n,m) ___ 

-°n0 5 -^nl 7***5 nm f 

with that precise ordering, and the interpolation nodes 

{(zi,J/i)|* = 0, ••• ,n;j = 0, ••• ,m} = 

{(£ 0 , 2 / 0 ), (^0,2/l), • • • , (^0,2/m), 

(^T, 2/o), (^1, 2/l), , (*^1, 2/m), , (•£«, 2/o), , (^n, 2/m)}, 

then the (n + l)(m + l) interpolation conditions F(xi,yj ) = fij can be written 
as a linear system 


dlc= /, 

where the coefficient matrix A is given by a Kronecker product 


with 


and 


= {(pf\xi)), i = 0,..., n; j — 0,... ,n, 
B v = ((^(2/i)), * = 0, • ■ ■, m- j = 0,... ,m, 


£ (^00 7 * * * 7 Q)m? ^lO? * * * 7 ^lm? 7 7 £nm) 7 


/ — (/0O7 * * * 7 /om, /1O7 * * * 7/1 m, , /nO, , fmn) 


The Kronecker product D ®E is defined by blocks as ( dkiE ), with D 


( dki )• 


For reasons which will be explained in Section 4 we will select as interpolation 
nodes (x i: yj) = (^|, ^) (i = 0, • • • , n; j = 0, • • • , m). In the general case 
we must avoid the value of x t for which the leading coefficient of p(t) in the 
power basis evaluates to 0, and the value yj for which the leading coefficient 
of q(t) in the power basis evaluates to 0. 
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In our example we have n — 3 and m = 4, and consequently B x will be the 
matrix 


B x 


/ 64 

48 

12 

_L\ 

125 

125 

125 

125 

27 

54 

36 

8 

125 

125 

125 

125 

8 

36 

54 

27 

125 

125 

125 

125 

. _L 

12 

48 

_6£ 

\ 125 

125 

125 

125/ 


and By will be the matrix 


/ 625 

125 

25 

5 

1 \ 

1296 

324 

216 

324 

1296 

16 

32 

8 

8 

1 

81 

81 

27 

81 

81 

1 

1 

3 

1 

1 

16 

4 

8 

4 

16 

1 

8 

8 

32 

16 

81 

81 

27 

81 

81 

1 

5 

25 

125 

625 

\ 1296 

324 

216 

324 

1296 ) 


As it is well known, since B x and B y are nonsingular matrices the Kronecker 
product B x ® By will also be nonsingular. 


As for the generation of the interpolation data, let us remark that they can 
be obtained without constructing the symbolic Bernstein-Bezout matrix BS. 
That is to say, we can obtain each interpolation datum by means of the eval¬ 
uation of p{t) and q(t) followed by the computation of the determinant of the 
corresponding numerical Bernstein-Bezout matrix B making use of the Bini- 
Gemignani algorithm which constructs (in 0(n 2 ) arithmetic operations) the 
Bernstein-Bezout matrix for the evaluated polynomials. 

In addition, we must divide the value of the determinant by —2 Xi + 1. 


An algorithm for solving linear systems with a Kronecker product coefficient 
matrix is derived in a self-contained way (in a more general setting) in [14]. 
For the case of the power basis considered in [12], taking into account that 
every linear system to be solved was a Vandermonde linear system, it was 
convenient to use the Bjbrck-Pereyra algorithm [3, 9] to solve those linear 
systems. For the Bernstein basis being used here, an appropriate algorithm 
which takes advantage of the special properties of the coefficient matrices B x 
and By will be presented in Section 4. 

In the general case, we must solve n + 1 linear systems with the same matrix 
By and m + 1 linear systems with the same matrix B x . 


7 








4 Total positivity of B x and B y 


From [4] we know that the Bernstein basis of the space of polynomials of 
degree less than or equal to n is a strictly totally positive basis on the open 
interval (0,1), which implies that all the collocation matrices 


M = (Pj n \ti)),i,j = 0 

with to < ti <...< t n in (0,1) are strictly totally positive, i.e. all their minors 
are strictly positive. In particular, due to our choice of the interpolation nodes 
the matrices B x and B y are strictly totally positive matrices. 

Making use of the results of [7, 8], we know that performing the complete 
Neville elimination on a strictly totally positive matrix A a bidiagonal factor¬ 
ization of its inverse A -1 can be obtained, that is to say, we have 


yr 1 = GiG 2 ... Gn^D-'Fn^F^ ...F,, 
where D is a diagonal matrix and F* and Gi are bidiagonal matrices. 

So, after having obtained that factorization (with a computational cost of 
0(n 3 ) arithmetic operations), all the systems Az = b with coefficient matrix 
A can be solved (with a cost of 0(n 2 ) arithmetic operations) by performing 
the product 


G\G 2 - ■ - G n -\D 1 F n _ 1 F n _ 2 ... Fib. 

An early application of these ideas to solve structured linear systems can be 
seen in [15], and a recent extension has been presented in [6]. 

A detailed error analysis of Neville elimination, which shows the advantages 
of this type of elimination for the class of totally positive matrices, has been 
carried out in [1], and related work for the case of Vandermonde linear systems 
can be seen in Chapter 22 of [10]. 

In our situation we must notice that the bidiagonal factorization can be done 
in exact arithmetic, and the results of the factorization can then be rounded if 
the subsequent computations must be carried out in finite precision arithmetic. 

After having obtained the bidiagonal factorization of the inverse of B y , the 
solution of the linear system B y z = b can be obtained in 0(n 2 ) arithmetic 
operations by computing the product 



G 1 G 2 ...G n - 1 D- 1 F n _ 1 F n - 2 ...F 1 b, 
and analogously for the linear systems with coefficient matrix B x [6]. 

In our example, the coefficients of the desired implicit equation in the tensor- 
product bivariate Bernstein basis (using the lexicographical ordering we are 
considering) are: 


(25264/27,66256/81,167852/243,45652/81,36137/81,15728/27,125312/243, 

320120/729,29164/81,69421/243,29440/81, 79024/243, 203228/729, 
18580/81,14761/81,2048/9,16640/81,14336/81,3940/27,9391/81). 


5 Computational complexity 


In this section we will briefly examine the computational complexity of our 
algorithm in terms of arithmetic operations. In view of the algorithm, we must 
solve n+1 systems of order m +1 with the same matrix B y and m +1 systems 
of order n + 1 with the same matrix B x . 


The factorization of the inverse of a matrix of order n by means of complete 
Neville elimination takes 0(n 3 ) operations, but that factorization is used for 
solving all the systems with the same matrix, so each of the remaining systems 
can be solved with 0(n 2 ) operations. 

For the sake of clarity in the comparison, we will consider here the case m = n. 
Then, the interpolation part of the algorithm has computational complexity 
0(n 3 ). Let us observe that in this situation, if we solve the linear system 
Ac = f of order (n + l) 2 by means of Gaussian elimination, without tak¬ 
ing into account the special structure of the matrix, we have computational 
complexity 0(n 6 ). Moreover, using the approach we are describing, there is 
no need of constructing the matrix A, which implies an additional saving in 
computational cost. 

Let us remark that, since the construction of the numerical Bernstein-Bezout 
matrix requires 0(n 2 ) arithmetic operations and the complexity of the compu¬ 
tation of each determinant is 0(n 3 ), the generation of the interpolation data 
has a computational complexity of 0(n 5 ). Therefore with our approach, which 
exploits the Kronecker product structure, the whole process has complexity 
0(n 5 ), while using Gaussian elimination it would be 0(n 6 ). 
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It is worth noting that the main cost of the process corresponds to the genera¬ 
tion of the interpolation data, and not to the computation of the coefficients of 
the interpolating polynomial. So, the main effort to reduce the computational 
cost must be focused on that stage. In this sense, an interesting issue would 
be to take advantage of the displacement structure of the Bernstein-Bezout 
matrices [2, 11] to develop an algorithm with complexity 0(n 2 ) for computing 
each determinant. 

Remark. Finally, let us observe that all the linear systems with matrix B y 
can be solved simultaneously, and the same can be said of the systems with 
matrix B x . Therefore the algorithm exhibits a high degree of intrinsic paral¬ 
lelism. This parallelism is also present in the computation of the interpolation 
data since we can compute simultaneously the determinants involved in this 
process. 
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